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Direct detection of vacuum fluctuations and analysis of sub-cycle quantum properties of the 
electric field are explored by a paraxial quantum theory of ultrafast electro-optic sampling. The 
feasibility of such experiments is demonstrated by realistic calculations adopting a thin ZnTe electro¬ 
optic crystal and stable few-femtosecond laser pulses. We show that nonlinear mixing of a short 
near-infrared probe pulse with multi-terahertz vacuum held modes leads to an increase of the sig¬ 
nal variance with respect to the shot noise level. The vacuum contribution increases significantly 
for appropriate length of the nonlinear crystal, short probe pulse durations, tight focusing, and 
sufhciently large number of photons per probe pulse. If the vacuum input is squeezed, the signal 
variance depends on the probe delay. Temporal positions with noise level below the pure vacuum 
may be traced with a sub-cycle accuracy. 

PACS numbers: 42.50.Ct, 42.50.Lc, 42.65.Re, 78.20.Jq 


Finite fluctuation amplitudes in the ground state of 
empty space represent the ultimate hallmark of the quan¬ 
tum nature of the electromagnetic radiation field. These 
vacuum fluctuations manifest themselves indirectly In a 
number of phenomena that are accessible to spectroscopy 
such as the spontaneous decay of excited atomic states as 
well as the Lamb shift [T] in atoms [5] and in quantum- 
mechanical electric circuits [3]. Access to the quantum 
aspects of electromagnetic radiation is provided by the 
analysis of photon correlation mi] or homodyning [i- 
111] measurements. However, these approaches require 
amplification of the quantum field under study to finite 
intensity and averaging of the information over multiple 
optical cycles. 

On the other side, precise determination of voltage or 
electric field amplitude as a function of time represents 
a fundamental task in science and engineering. Opti¬ 
cal techniques have to be applied when detecting electric 
fields oscillating in the terahertz (THz) range and above. 
Those approaches involve probing with ultrashort laser 
pulses of a temporal duration on the order of half an os¬ 
cillation period at the highest frequencies under study. 
Far-infrared electric transients [n na can be charac¬ 
terized by photoconductive switching [T3]. Electro-optic 
sampling in free space [Mil! allows field-resolved detec¬ 
tion at high sensitivity in the entire far- and mid-infrared 
spectral range [B ini. Direct studies of the complex¬ 
valued susceptibilities of materials and the elementary 
dynamics in condensed matter may be performed with 
these methods [BlU]- The time integral of near-infrared 
to visible electric-field wave packets is accessible with at- 
tosecond streaking [22]. So far, all those techniques were 
restricted to the classical field amplitude. 

In this Letter, we demonstrate theoretically that the 
quantum properties of light may be accessed directly in 
the time domain, i.e. with sub-cycle temporal resolution. 
Our considerations are based on the realistic example 
of electro-optic detection with zincblende-type materials 



FIG. 1. (color online) Electro-optic sampling setup and geom¬ 
etry. (a) The incoming near-infrared (NIR) probe and multi- 
THz signal fields mix in the electro-optic crystal (EOX). The 
NIR (blue) spatial mode amplitude is depicted by the con¬ 
tour plot whereas a sampled THz (red) spatial mode is in¬ 
dicated by wave fronts. Bottom left corner: time profiles of 
the NIR intensity envelope /nir(^) a representative multi- 
THz vacuum field E,^-^^JA). After collimation by a lens (L), 
the modified NIR field is analyzed using a quarter-wave plate 
(^)^, a Wollaston prism (WP) and balanced detectors (Ds,Dz) 
measuring the difference in the photon flux for the split com¬ 
ponents. (b) Spatial directions determining the electro-optic 
effect in the zincblende-type EOX and the following ellipsom- 
etry analysis. 


[23] . Even vacuum fluctuations may be sampled without 
amplification by broadband probing of electric field am¬ 
plitudes in the multi-THz region with few-femtosecond 
laser pulses of moderate energy content. 

We consider the geometry of electro-optic sampling 
shown in Eig.[^ An ultrashort near-infrared (NIR) wave 
packet with electric field Ep propagates along the [110] 
axis of an electro-optic crystal (EOX) |23l [24] . Its wave 
vector is perpendicular to the z-axis of the EOX. 
We select Ep || In this configuration, the second- 

order nonlinear mixing of Ep(t) with an incident THz 
field E.j,jj^(t) induces nonlinear polarization in the EOX 
plane with the components (for details, see Ref. [35] ) 

=-eodE,^Ut)E^it), Pi^Ht)=0. ( 1 ) 
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eo is the vacuum permittivity. The coupling constant 
d = can be determined from the electro-optic 

coefficient and refractive index (RI) n at the cen¬ 
tral frequency Wc of Ep [5SH1H]. In general, both fields 
= Ethz.s and in Eq. 0 are quantized, whereas 
Ep = Ep^z = {Ep,z) denotes the classical part of the 
probe held. We neglect the effect of quantum mechan¬ 
ical huctuations of the probe held on assuming a 

sufficiently large Ep. 

The nonlinear polarization generated by the wave 
mixing in the EOX represents a source in the inhomoge¬ 
neous wave equation describing propagation of the elec¬ 
tric held E in the EOX. The helds F = E, P*^^^ prop¬ 
agating in the forward direction (see Fig. can be 

decomposed as F(r,t)=J'((^dw F(r; where 

= wn^/cQ. Co and are the velocity of light and the 
frequency-dependent RI of the EOX, respectively. Using 
the paraxial approximation EiEn], the inhomogeneous 
wave equation reads 

E(r;cc) = -^p(2)(r;a;), (2) 

eoco 

where and ^ -f From Eq. Q 

we obtain Pi^^(r;w) = —eod f^di2 E.^^^ (r; fl)Ep(r; uj — 

^ electric held of the probe beam 
provides a solution of the homogeneous part of Eq. 0 
which can be decomposed into Laguerre-Gaussian (LG) 
modes isiisa (see Ref. [l5]). We adopt a probe pulse 
train with a fundamental Gaussian transverse mode of 
amplitude Q;p(a;): 

Ep(r; w) = ap(a;)LG„„ (r^, r,,; . (3) 

A length I of the EOX much shorter than the Rayleigh 
range of a beam at the relevant THz frequencies U with 
waist size wq is assumed, i.e. I A:nWg/2. 

The EOX is located at the beam waist, = 0, and 
has anti-rehection coating on its surfaces. Denoting 
F(r^;a;) = F(r^,r|| =0;w) we hnd that at the exit from 
the EOX, r|| = 1/2, the total electric held in the (110) 
plane is given by 

E'{y) = Ep{y)ez + {y)es + SE\y) , (4) 

where ^ = {rj^;tLi}. 5E'(y) = Ep/y) — Ep{y)ez denotes 
the vacuum held contribution at the probe frequency w 
in the vacuum picture [33]. The correction to the probe 
held generated in the EOX is evaluated as 


A 


2ika 


dr„ 


is the RI at U, whereas n and rig are the RI and the 
group RI codki^j/duj at w = Wc, respectively. Going be¬ 
yond Ref. |T3 where an expression similar to Eq. ([^ was 
derived for the case of plane waves in order to establish 
a classical theory of electro-optic sampling, Eqs. 0 and 
([^ include the transverse spatial dependence of the helds, 
the quantized form of the signal as well as the contribu¬ 
tion of quantum huctuations at the probe frequencies. 
These points are crucial for our further analysis. 

From Eq. 0: we see that the nonlinear mixing of the 
probe and THz components generates a new held prop¬ 
agating in the same direction and polarized perpendic¬ 
ular to the probe. For the analysis of the polarization 
state of the modihed probe, we consider the held com¬ 
ponents in the coordinate frame ea = (e^, eg)/y/2 ro¬ 

tated by 45° with respect to the frame [Fig. [^b)], 
&a{y) = Ep{y)[l±i/>{y)]/V2 + 5&a{y) . Here ^{y) = 

b b 

iE^^'> {y)/Ep{y) must be small for the frequency range 
of the probe. 

The ellipsometry setup used in typical experiments is 
explained in Fig. Sa). We consider its ehects at the exit 
surface of the EOX. This simplihcation is justihed when 
all probe photons are detected without spatial hltering. 
The hrst step of the analysis consists in describing the 
action of the quarter-wave plate with axes oriented along 
Ba and Gh such that it phase-shifts the a-component of 
the held by n/2: E'/(y) = tE/iy), &/{y) = E',{y). 
The Wollaston prism splits the electric held into its z- 
and s-components: 

^{y) = [I T ky)] + s&ziy ). (6) 

Finally, the photon numbers in both held components are 
detected and subtracted. The photon number operator 
for the polarization a = z,s reads |34| 

= G^“dw^ J dV,E"t(r,;a;)E"(r,;a;) , (7) 

where C = dTrcgiieo, the dagger denotes Hermitian 
conjugation and the spatial integral covers the entire 
transverse prohle of the probe beam. The frequency- 
dependent quantum efficiency of the photodetector 
ri{ui) ~ I over the detected frequency range but vanishes 
quickly for w —>■ 0. 

Inserting Eq. Q into Eq. 0 and neglecting the 
second-order terms in dE" as well as the mixed terms 
depending linearly both on (5E" and on (contained 
in (j)) [5^, we obtain for the total detected quantum sig¬ 
nal 


/ OO 

-OO 


where the factor = “^^^(^sinc 


-%) 


de¬ 


termines phase matching. Here sinc(a;) = sin(a:)/a;, 


S = //'s ~ — Seo + ^sn , 

where the electro-optic signal (EOS) Seo is 
r]{uj) 


( 8 ) 


5eo = C / duj^ |Ep(>;)|^ c/iy) + H.c. (9) 
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and the shot 

= Cfd\f^ 


noise (SN) 


0 nu. 


contribution 5sn 
E;(y)SEUy)+ H.C. 


reads 

Here 


H.c. denotes the Hermitian conjugate and SE'^(y) = 
e^'^/'^lSE'J(y) + iSE"(y)]/v^ is the circular component 
of the probe field vacuum contribution |36j . Summing up 
the signals from both detectors, we obtain the expecta¬ 
tion value of the number of detected photons per probe 
pulse N = (A^s+K) = ^^Y^rdu;^lap(u;)l^. 

Using Eqs. (§ and © in Eq. we obtain 


~ _ dlNujp 

^eo — 

Con j 


-5oo(u) 


poo 

/ dUE^ 

' —oo 


K;U)i?(U). (10) 


£'oo(u) = LGj,„(rj^,r|| =0;fc^) = w/exp(-r2/u;2) 

is a normalized Gaussian independent of uj and Wp = 
/“dw 77(a;)|ap(a;)|Y/“da;^|ap(a;)|2 is the aver¬ 
age detected frequency. We have introduced the re¬ 
sponse function i?(U) = sine — Ug) f{fl) with 

the normalized Hermitian spectral autocorrelation func¬ 
tion /(H) = [/*(H) -H /_(H)]/2, where /^(H) = 

ir al{uj)ap{u} ± H)//“ duj r]{uj)\ap{uj)\‘^ . 
Within the paraxial quantization [i 


in 


Eq. (10) is given by 





hn 


dTreoCoUs 


-as,i.p(f^)5' (U) (11) 


for H > 0, < 0) = El^^{r^;-V,). Here, 

ns,i,p(U) annihilates a photon with frequency H, or¬ 
bital quantum numbers l,p and polarization eg. We 
have introduced the transverse mode functions (r^) = 
LG,j^(r^,r|| = 0;fcn). In contrast to the probe beam, the 
waist size Wg characterizing these mode functions is a 
free parameter of the expansion ( |11[ ). Inserting Eq. 0 
into Eq. (10) and selecting Wg = wg/'/2, we can per¬ 
form the spatial integration using f d^r^g^^ (r^)g'^ (r^) = 
^wo ^‘.o^p,o- Then we obtain from Eq. (10) 


Seo = -iVS 



[ag,o.o(ll)i?(ll)-H.c.], (12) 


where B = (d^l^N^ujpfi)/(47r^eoCgn^Wg). 

As an input, we now consider a THz quantum field 
with no coherent (classical) contribution: (E^hz) = Oj 
e.g., a bare multi-THz vacuum. Then (S) = 0 since 
(Ssn) = 0 and ^ in Eq. © depends linearly on E^hzj 
thus also (Seo) = 0. However, the variance of the sig¬ 
nal does not vanish. If the range of detected THz fre¬ 
quencies, determined by i?(H), does not overlap with the 
frequency content of the probe beam, the signal variance 
(5^) —(5)^ = (S^) can be written as (S^) = (Sg^) + (S^n)- 
Calculating the SN contribution using the paraxial quan¬ 
tization [32], we obtain the expected result (Sg^} = N. 



(color online) (a) Calculated integrand function 

(b) Double-logarithmic 


FIG. 2. 

U (n/n^)|i?(U)P entering Eq. (|13 1. 


plot of the ratio AS/N in dependence on N. Black dotted 
(red dashed) line shows the bare SN (multi-THz vacuum) con¬ 
tribution. (c) Increase of (A5 — A5sn)/A5an with N. Param¬ 
eters are defined in the main text. 


Evaluating (5eo) multi-THz vacuum yields 


Cq 


lujp Y fi/“dHH(n/nJ|i?(H)|2 


dlT^egCgnWg 


(13) 


where we have used (ag^o,o(f^)aI o o(^0) = d(H —H'), 
whereas the expectation values of other possible 
quadratic combinations of ^ and dg^i^p vanish. Note 
that the second and third factors on the right-hand side 
of Eq. (13) have the dimensions (m/V)^ and (V/m)^, re¬ 


spectively. The latter can be interpreted as the square 
of the effective multi-THz rms (root mean square) vac¬ 
uum electric field filtered by the response function. The 
former, oc , determines how effectively this field is sam¬ 
pled for a fixed N. 

To illustrate the results, we assume the following real¬ 
istic specifications of the sampling few-femtosecond NIR 
laser pulse: center frequency 255 THz, spectral band¬ 
width 150 THz with rectangular spectral shape and flat 
phase, leading to Wp = 247 THz, and waist size Wg = 
3 p,m [37]. We consider a Z = 7 /rm thick ZnTe EOX 
with = 4 pm/V [351 I3S|, n = 2.76, rig = 2.9, and 
varying only slightly (from 2.55 to 2.59) for relevant 
THz frequencies |25| . The resulting integrand function 
entering Eq. (13) is shown in Fig. [^a) (for details, see 
Ref. [IS]). Diffraction effects are taken into account by 
excluding wavelengths A with X/{2n^) > wg. 

Based on this input, we calculate the dependence of 
the rms value of the signal AS = (5^)^^^ on the aver¬ 
age number N of photons in the sampling NIR pulse, as 
shown in Fig. [^b) on a double-logarithmic scale. Above 
a certain N, the EOS contribution of the multi-THz vac¬ 
uum changes the typical SN scaling. The relative increase 
of the rms value of the signal with respect to the SN level, 
{AS — ASsn)/ASsn, is depicted in Fig. [^c) for moder¬ 
ate N and with linear scaling. For even higher N, the 
vacuum contribution starts to dominate so that the de¬ 
pendence saturates to the constant EOS level [Fig. [^b)]. 
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FIG. 3. (color online) (a) Frequency dependence of the 
squeeze factor r. Squeezing correlates fl and 2flc — 
modes (as indicated by arrows), (b) Error contour in the 
complex-amplitude plane for PV (gray circle) and SV with 
9 = 0/8 = n (red/green ellipse with reduced uncertainty in the 
phase/amplitude quadrature YjX). (c) Normalized (with re¬ 
spect to the constant PV level, solid black line) EOS variance 
in dependence on the time delay r of the probe NIR pulse for 
SV with 9 = 0/6 = 11 (solid red line / dashed green line). 


Subtracting the SN contribution from the total signal 
variance, the bare EOS variance induced by the sampled 
quantum field can be obtained and analyzed. 

To elaborate on this point, we apply our the¬ 
ory to a multi-THz vacuum which is squeezed in 
an interval around a central frequency Oc. The 
corresponding state of light is obtained by ap¬ 
plying the continuum squeezing operator |40lH^ 

exp \ /g dfl Y.aip - Vt)aa,i,p{^) - H.c.) 

to the multi-THz pure vacuum (PV) state considered 
above. Here the frequency-dependent squeezing param¬ 
eter ^(H) satisfies the condition ^(H) = ^(2Hc — H). 
We assume that all spatial and polarization modes 
are squeezed equally. In this case, the EOS can be 


obtained from Eq. (12) applying the transformation 
hs,o,o(f^) 


as,o,o{^) coshr^ — aj q o(20c — sinbr^^ 

and working in the vacuum picture. The ex¬ 
pectation value of the signal remains zero. Evaluation 
of the EOS variance for the squeezed vacuum (SV), 
(5eo)sv(T), is analogous to the PV case. However, 
the SV EOS variance depends on the time delay r 
of the NIR probe pulse leading to the transformation 
i?(H) —>■ ii(H)e“®^'^ of the response function, a fact that 
was unimportant for handling the PV [cf. Eq. (HU)]- 
For a probe pulse symmetric with respect to t = t, i.e. 
Ep{t — t) = Ep{T — t), we find R{Q) = , 

where Ro{^) is real-valued. 

For our illustration we assume constant squeezing with 
^(H) = ^ = re*® in the frequency range with 

He = (f2i -b H 2 )/ 2 , where r = |^| is the squeeze factor 
[inilll] and 9 = Arg(^) is the squeezing phase [41]. No 
squeezing occurs outside this range. In particular, we 
use ^/(27r) = 40 THz, H 2 — Hi = He and sinhr = 2 [see 
Fig. ^a)]. Generalized quadrature operators [ID] X\ = 


^ 2 ( 02 - 01 ) [o^,o,o(^)e *^-bH.c.], with X = Xo 

and V = Xj^/ 2 , normalized so that [V,?] = i are in¬ 
troduced. The error contours for PV as well as for SV 
as described above and two different squeezing phases, 
9 = 0 and 0 = tt, are featured in Fig. I^b). The depen¬ 


dence of the normalized EOS variance {S^o)sv{t)/{S^ o), 
where (iSgo) is given by Eq. (13), on the time delay t is 
shown in Fig.|^c) for the same states as in Fig. [^b) and 
sampling parameters used for the PV case. For specific 
delay times, the EOS variance of the multi-THz SV is 
by 64% lower than the unsqueezed value of the PV state 
(for details, see Ref. pSj). 


We emphasize the cardinal difference between our find¬ 
ings and similar-looking results obtained in the context 
of homodyning |45ll46j . In homodyning experiments, the 
signal is determined by the temporal overlap integral of 
the complex amplitudes of the electric fields of an input 
state and a local oscillator, i.e. the information is essen¬ 
tially averaged over multiple oscillation cycles of light. A 
restricted frequency bandwidth has to be assumed to jus¬ 
tify the slowly varying amplitude approximation under¬ 
lying this approach. In contrast, electro-optic sampling 
provides a true sub-cycle resolution of the probed multi- 
THz electric field. Moreover, registration of photons is 
transferred into the NIR, circumventing the lack of effi¬ 
cient single-photon detectors in the multi-THz frequency 
range. Most importantly, the multi-THz quantum field 
may be studied without the necessity to reduce or am¬ 
plify its photon content - even if it remains in its ground 
state. 


In conclusion, we theoretically clarify the contribution 
of the quantum fluctuations of the multi-THz vacuum 
electric field to the signal in ultrabroadband electro-optic 
sampling by differentiating it from the trivial shot noise 
of the high-frequency gating pulse. The crucial aspects 
are a strong localization of the sampling beam in space 
and time as it passes the nonlinear crystal, a large second- 
order nonlinear coefficient and proper phase matching 
that might be further optimized selecting an even more 
appropriate material than the thin piece of ZnTe we have 
considered as an example. For a multi-THz squeezed vac¬ 
uum, the possibility to trace the oscillations of the EOS 
variance with the time delay of the probe pulse is pre¬ 
dicted. Positions occur where the noise remains signifi¬ 
cantly below the level of unsqueezed vacuum. The same 
formalism can be applied for the analysis of more complex 
quantum fields in a time-resolved and non-destructive 
manner. Experimental implementation of these ideas 
might open up a new chapter of quantum optics oper¬ 
ating predominantly in the time domain and with access 
to sub-cycle information on the quantum state of electro¬ 
magnetic radiation. 
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Supplemental Material 


Paraxial Theory of Direct Electro-Optic Sampling of the Quantum Vacuum 

A.S. Moskalenko, C. Riek, D.V. Seletskiy, G. Burkard, and A. Leitenstorfer 


1. GEOMETRY OE NONLINEAR MIXING 


For the nonlinear mixing in the BOX, we select the orientation of Ep parallel to the z-axis. This choice of the 
polarization direction of the probe field ensures that the maximum signal is detected in the electro-optic detection 
scheme for a copropagating classical THz electric field polarized perpendicular to the probe electric field [ST]. In the 
experiment, adjustment is achieved by rotation of the BOX around the [110] axis for fixed, mutually perpendicular 
polarization directions of the probe and detected electric fields. Also, only one of two possible polarization modes of 
the detected field [the one perpendicular to the probe field, i.e. oriented parallel to the unit vector in Big. [^b)] 
contributes to the signal in this geometry. There is no THz field generated by optical rectification of the probe for 
this orientation of the BOX. 

The second-order nonlinear mixing of the probe field Ep(t) with the detected THz field then induces 

nonlinear polarization in the BOX with the following components [S2HS4) : 

Px = ^^odxyzE^^^^y{t)Ep^z{t) (SI) 


and similarly for the y-component with the interchange of indices a: -O- y in Bq. (SI). Here eg is the vacuum permittivity. 
For the zincblende-type FOX we adopt as an example, the tensor components d^yz = X^yzl‘^ and dyxz = X'yxzl‘i‘ are 
both equal to the same constant denoted by d^^ [S2HS4] . This coefficient is related to the constant used for the 
description of the Pockels effect as d^^ = —n^r^j/4. n is the refractive index at the central frequency of the probe 
electric field. For the following discussion, it is convenient to introduce d = Writing Bq. (SI) as 

an instantaneous relation in the time domain we assume that the frequencies H of the THz field are lower than the 
frequencies oj of the probe field and that the second-order nonlinear coefficient can be considered constant in the 
frequency range determined by the spectral width of the probe field. The frequency dependence of the nonlinear 
coefficient can be easily included writing the corresponding equations in the frequency domain, similar to Ref. |S5j 
but taking care of the particular geometry. However, the discussion of the geometrical issues is more concise in the 
time domain whereas the effect of the frequency dependence of the nonlinear coefhcient is finally not significant in 
our case. 

The nonlinear polarization induced in the (110) plane is given by = '^{Py^'^ ~ 


= -^(Py^’ — Px‘‘’), using the unit vectors eg and as a basis in this plane. Taking Bq. (SI I, expressing also 
the components of the quantized THz field in this basis and neglecting the effect of quantum mechanical fluctuations 
of the probe beam on the induced nonlinear polarization, i.e. assuming a sufficiently strong probe field, we arrive at 
Bq. (0. Inclusion of the vacuum contribution for the probe beam at this place would mean taking into account mixed 
second-order corrections linearly dependent on both the vacuum fluctuations of the probe field and on the probed 
THz field. In our present consideration, we neglect such terms since they do not lead to significant effects. 
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= 0 and 


2. LAGUERRE-GAUSSIAN MODES AND PARAXIAL ELECTROMAGNETIG FIELD QUANTIZATION 


In electro-optic sampling, propagation of the NIR probe beam through the BOX can be well described within the 
paraxial approximation. The same approximation can be naturally used to describe the sampled multi-THz quantum 
fields. The corresponding expression for a quantized electric field within the paraxial approximation was derived in 
Ref. [56]. In free space, with the propagation axis selected as shown in Fig. it reads 


E(r,t) = -f XI 
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where da^i^p{k) denotes the annihilation operator for a photon with absolute value of the wave vector fc, frequency 
Q = cofc, orbital quantum numbers Z,p, and polarization direction eo,. The spatial mode functions are given by the 
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Laguerre-Gaussian (LG) modes LGip(r^, r,,; k) = LGip(rj^, ip, r,.; k) which can be written as 


LGip{r^,ip,r ;k) = 


2p\ 


'v^r 


W 


tt{\ 1\ +py. wir ) ) w{r ) 




2r^ 


'(g). 


exp 


2r^ 


'(g) 


■ Up 


kr^ 


'27^(r„) 


• (S3) 


Here rc(r||) = ico+ ’’^/^r(^) is the transverse mode radius at the longitudinal position with wq being the 
waist size of the probe beam (mode radius at = 0) and Ir{^) = fcrcg/2 denoting the Rayleigh range of the 
beam at given k. = r|| l + Z^(n)/r^ is the phase-front radius, d?G(?'||) = —{‘^P+ |^| + 1) arctan(r||/wg) is 

the Gouy phase and Lp^{x) are the associated Laguerre polynomials [S8]. The LG modes are normalized such that 
/g’^df/)/g°°dr^r^ LG*p{rj^,4>,r^^-k)LGi'p'{r^,4>,r^^',k) = SwSpp' (for any k and r||), where Sij denotes the Kronecker 
delta. The annihilation and creation operators satisfy the continuum commutation relations [aaj^p{k),aapii^pi{ky] = 
[®a i p(^)’ ^a' i> p' (^0] = 0 [^a,i,p{k),dl^, p pi{k')] = 5aa'dw5ppi5{k — k'). Expressing the creation and annihilation 

operators as functions of frequency, whereby they satisfy [aQ,,/,p(n),p'(fl')] = i p{^)tdli, p p,{k')] = 0 and 
[aQ,i,p(f2), p ,(H')] = 5aa'Sw5pp'5{il. - il'), Eq. (p^ transforms into 


E(r,t) = -iJ2 
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(S4) 


By writing k^ we have stressed that we consider k = k^ = fl/cg as a function of H in this expression. Note that 
the factor of 167r^ in the denominator under the square root in Eq. (20) of Ref. [S6] needs to be replaced by 47r. 
This fact is confirmed by deriving the expression for the total energy operator of the electro-magnetic field and 
has been considered in Eq. ([S2l). From Eq. ([S4|, the total energy operator E is obtained in its correct form as 
£ = /“dG fif2E„,.p«I.z.p(f^)Sa.z,p(f^). 


In media with refractive index n^, the factor under the square root in Eq. (S4) should be additionally divided by 
[Ml pp. 391-392]. This measure again ensures a correct expression for the total energy operator of the field so 
that Eq. (S41 takes the form 


E(r,t) = -iJ2 
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When we consider a thin EOX located at the beam waist (r|| = 0), we use LGip(rj^, Th ; fc) 
gip{r^), which are given by 


(S5) 
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The transverse modes gip{rj^) = gip(r^, </)) are independent of k. They are normalized such that 


^277 pOO 

/ d()) / dr^r^ glj^(r^,4>)gpp,{r^,4>) = Sip6py . 

Jo Jo 


The fundamental (lowest-order) mode is Gaussian-shaped and given by 


9ooiT^±) = 


2 1 
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When the EOX is located at 0, the same approximation applies and one can proceed similarly, just changing 
Wo —> rc(?'||) in the expressions for gip{r^). 


3. REFRACTIVE INDICES AND RESPONSE FUNCTION 


In the paper we used a sampling few-femtosecond NIR laser pulse of the following specifications m-- center 
frequency uJc/{ 2 tt) = 255 THz, spectral bandwidth ALLi/(27r) = 150 THz with rectangular spectral shape and flat 
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FIG. SI. Temporal profile of the intensity (red line) and its envelope (blue line) of the NIR probe pulse used in our calculations. 
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FIG. S2. Dispersion of the refractive index in the THz (solid red line) and NIR (dashed blue line) range. Inset shows the 
dispersion with a higher frequency resolution in the range of small THz frequencies. 


phase. The corresponding temporal profile of the NIR probe intensity is shown in Fig. |S1[ For such a pulse we get 
Wp = 247 THz, where Wp is defined in the text after Eq. (10). Notice that Wp « Wc. However, there is a small difference 
in these quantities due to different averaging used in their definitions. This difference is of minor importance for our 


consideration. The normalized Hermitian spectral autocorrelation function /(H), as defined in the text after Eq. (101, 
can be found as 




(S8) 


where H(x) denotes the Heaviside step function. For our example with a rectangular probe spectrum, /(H) takes the 
shape of an isosceles triangle with the vertex at H = 0. In order to determine the response function i?(H), we have to 

multiply /(H) by the phase-matching function sine — rig) . The latter requires knowledge about the refractive 

index n in the THz range and group refractive index rig at the (NIR) central frequency Wc of the probe pulse. 


Refractive index properties of a ZnTe crystal in the NIR frequency range are modelled by a Sellmeier formula [Sll] 

nl=A+[B\^l{X^-c^)l (S9) 

with A = 27rco/H, A = 4.27, B = 3.01, and (? = 0.142 |S12j . The corresponding frequency dependence is shown on 


the right side of Fig. S2 For the refractive index in the THz frequency range, we use the parametrization from 
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frequency Q.l{2%) (THz) 

FIG. S3. Calculated response function R(yi) without the low-frequency cntoff (solid black line) and with the low-freqnency 
cutoff (dashed red line). 


Ref. [SH]: 


= Re 
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(SIO) 


— 206 cm 


7 = 3.01 cm“^, and Coo = 6.7. The corresponding frequency dependence 
From these models we obtain the following values of the refractive index and 
the group refractive index at oJc =255 THz: n = 2.76 and rig = 2.9. It is important that these indices are almost 
constant in the neighborhood of uJc- Using the calculated and rig in the phase-matching function with I = 7 /rm, 


with = 177 cm 
is shown on the left side of Fig. S2 


the response function R{VL) depicted in Fig. S3 results. As discussed in the text of the paper, we introduce a low- 


frequency cutoff excluding wavelengths A with X/{2n^) > rco, in order to take into account diffraction losses. The 


modified response function is also shown in Fig. S3 The resulting integrand function entering the integral in Eq. (13) 
is found in Fig. [^a) of the paper. Note that without introducing the low-frequency cutoff we would get just a small 
increase of approximately 21% for the EOS variance calculated from Eq. (131. 


4. EOS VARIANCE FOR THE SQUEEZED MULTI-THZ VACUUM 


For the EOS variance {S'^o)s 
Eq. (13), we obtain 


,(t) normalized with respect to the EOS variance of the pure vacuum given by 


0l)sAr)/{Sl) = l+2aM+2b^M{M + 1) cos(0 - 


(Sll) 


Here M = sinhr, a = la/1, and b = 1^/1, where 


/ = 
Ia = 

h = 


dO g^{n), 


dO 


dO p(Q)p(2r2(. — Q), 


(512) 

(513) 

(514) 


with p(Q) = The real coefficients a and b generally satisfy the Cauchy-Bunyakovsky-Schwarz inequality 

b < a. Obviously, also a < 1 is valid. The dependence of (5go)sv('r)/(5eo) oa the time delay r is shown in Fig. |^c) 
for the same states as in Fig. j^b) and sampling parameters used for the pure vacuum case. We clearly see that 
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for certain delay times the EOS variance of the squeezed multi-THz vacuum can beat the uncertainty limit set by 
the pure vacuum state. Here, for M = sinhr = 2, its minimum value constitutes « 36% of the pure vacuum level. 
For the selected parameters, a slightly stronger suppression of the quantum noise down to 34% can be achieved by 
increasing M to « 5.8, whereas the maximal noise is more than doubled. These values are generally determined by 
the coefficients a and b, introduced above, for which we have 6 < a < 1 in the considered case. Here, due to the 
limitation set by this inequality, a complete noise suppression in the EOS variance is impossible for any time delay. 
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